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TECHNICAL  MEMORANDUM  87/209 


CanadS 


Abstract 

Methods  for  calculating  the  moments  of  an  arbitrary  polygon  and  for  determining  whether 
a  point  lies  within  a  polygon  are  derived  and  discussed.  The  methods  are  efficient,  robust, 
concise,  and  easily  programmed  in  any  computer  language.  A  FORTRAN  77  subroutine  which 
calculates  the  first  three  moments  of  an  arbitrary  polygon  is  also  included,  as  is  a  subroutine 
which  determines  whether  a  point  lies  in  an  arbitrary  polygon. 


On  traits  das  ndtbodas  panssttant  de  calculer  las  nonants  d'un  polygons 
arbitrairs  at  da  ddtarniner  si  un  point  sat  situ*  *  l'intdrisur  d'un 
polygons.  Las  ndthodes  sont  efficsces,  solidss,  concises  at  facilaa  * 
prograaamr  dans  n* inporta  quel  langage  infornatiqua.  On  prdsanta  auasi  un 
sous-progranne  on  FORTSAI  77  qui  calcule  las  trois  praaders  nonants  d'un 
polygons  arbitrairs,  da  nlas  qu'un  sous -progresses  qui  ddtarnins  si  un  point 
sat  aitu*  *  l'intdrisur  d'un  polygons  arbitrairs. 
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Notation 


The  distance  of  the  origin  to  the  £**  side  of  the  polygon. 
The  n***  moment  tensor  for  the  polygon. 


rf’.jJ1’: 

Jfi’.jg’.lg’: 


Components  of  the  first  moment  of  the  polygon. 


Components  of  the  second  moment  of  the  polygon. 


The  largest  integer  which  does  not  exceed  n. 


Outward  pointing  unit  normal. 


The  number  of  sides  of  the  polygon. 


7): 

agn(*): 


Tensor  permutation  function  defined  in  Section  2.1. 


Function  defined  in  equation  (S.l). 


Position  vector. 


n  times 

The  tensor  dyadic  Soot .  .x. 


Ardength  around  the  polygon  perimeter. 


Variable  used  to  parameterise  a  side  of  the  polygon. 


Coordinates. 


Unit  vectors  in  the  coordinate  directions. 


Xk,Vk-  Coordinates  of  the  kth  vertex  of  a  polygon. 


Dirac  delta  function. 


Bold  face  characters  are  reserved  for  vectors  and  tensors. 
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1  Introduction 


A  problem  that  arises  in  many  engineering  applications  is  to  find  some  moment  (e.g  the 
area  or  the  centroid)  of  an  arbitrary  polygon.  At  DREA  the  problem  has  arisen  in  the  context 
of  the  calculation  of  potential  flows  via  panel  methods[lj.  At  large  distances,  the  potential  due 
to  a  panel  is  determined  most  efficiently  by  a  multipole  expansion  in  which  the  moments  of 
the  polygonal  panel  appear.  In  this  memorandum  an  efficient,  robust,  and  easily  programmed 
method  for  calculating  the  moments  of  a  polygon  is  derived.  Although  it  seems  likely  that 
the  method  has  been  derived  previously,  it  does  not  appear  to  be  widely  known.  Indeed,  the 
problem  was  considered  suitable  (though  it  was  not  discussed)  for  a  seminar  at  the  Mathematics 
Dept,  of  Dalhousie  University  in  which  problems  of  unknown  solution  were  to  be  tackled[2]. 
Further  evidence  comes  from  the  potential  flow  program  EN967[3]  which  calculates  the  moments 
of  quadrilateral  panels  using  a  method  which  is  less  efficient,  less  robust,  and  less  general  than 
the  one  presented. 

A  simple  and  efficient  solution  to  a  related  problem  is  also  derived  in  this  memorandum: 
how  does  one  determine  whether  a  given  point  lies  inside  or  outside  an  arbitrary  polygon?  It, 
too,  has  arisen  in  the  context  of  potential  flow  panel  methods;  it  is  sometimes  necessary  to 
know  whether  a  certain  point  lies  within  a  panel.  This  problem  has  also  arisen  in  modelling  of 
a  "cycle  of  perception”  for  a  computational  vision  problem  by  the  Computer  Aided  Detection 
Group  at  DREA[4], 

Appendix  A  contains  a  FORTRAN  subroutine  which  calculates  the  first  three  moments  of 
an  arbitrary  polygon  using  the  method  discussed  in  this  memorandum.  Appendix  B  contains 
a  FORTRAN  subroutine  which  determines  whether  a  point  lies  inside  or  outside  an  arbitrary 
polygon. 


2  Calculation  of  the  Moments  of  Polygons 


In  this  section  the  method  of  calculation  of  the  moments  of  an  arbitrary  polygon  is  derived 
and  discussed. 


2.1  Analytical  Formulae  for  the  Moments 


Let  (x,  y)  be  a  coordinate  system  with  unit  vectors  £  and  y  along  its  axes.  Bold  face 
characters  will  be  used  to  denote  vectors  and  tensors.  Thus, 


_ _  A  ,  A 

x  =  xz  +  yy 


The  notation  xn  will  be  used  to  denote  the  tensor  dyadic 


n  times 


The  problem  to  be  solved  may  be  stated  as  follows: 


Problem  1:  Ifxk,k  =  1  are  the  vertieea  of  a  polygon  in  order  as  one  proceeds  around 

its  perimeter  counterclockwise,  calculate  the  nth  moment  of  the  polygon, 


l(")  =  [  xndxdy 

Ja 


where  the  notation  □  denotes  integration  over  the  surface  of  the  polygon. 


l(°)  is  the  area  of  the  polygon  and  iW  is  its  centroid  times  its  area.  If  the  polygon  has  uniform 
density,  the  second  order  tensor  is  proportional  to  its  moment  of  inertia. 

The  essence  of  the  method  is  to  express  the  xn  as  the  divergence  of  a  tensor,  so  that  the 
divergence  theorem  can  be  used  to  express  the  moment  as  a  line  integral  around  the  perimeter 
of  the  polygon.  The  contribution  to  the  line  integral  from  each  side  is  calculated  easily.  By 
using  tensor  notation,  one  can  obtain  a  single  expression  for  any  moment  of  the  polygon. 

First,  note  that 


V  •  xn  = 


d[xxn~l)  djyx"-1) 
dx  dy 


=  2xn_1 


dx”-1 


=  2xn  1  +  x  •  Vx" 
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Now,  since  x  •  Vx  =  x,  one  has  x  *  Vxn  =  nxn_1,  whence  from  equation  (2.4) 

V.x"  =  (n  +  l)x"-1  (2.5) 

Therefore,  using  equation  (2.5),  the  divergence  theorem,  and  the  definition  of  equation  (2.3), 
one  obtains 

I(n)  =  f  V  •  xn+1dxdy  =  [  h  •  xn+1d«  (2.6) 

n  +  2  Ju  n  +  2/sD 

where  denotes  the  perimeter  of  the  polygon,  n  is  an  outward  pointing  unit  normal,  and  da 
is  an  increment  of  arclength. 

The  kth  side  of  the  panel  may  be  parameterized  by  x  =  [(x*+i  +  x*)  +  t(x*+i  -  x*)]/2, 
t  €  [—1,1].  The  increment  of  arclength  is  then  da  =  |x*+i  -  x*|dt/2.  The  outward  pointing 
normal  is  parallel  to  (x*+i  —  x*)  x  £  where  £  is  a  unit  vector  perpendicular  to  the  plane 
of  the  polygon  and  such  that  £,  y,  and  £  define  a  right  handed  coordinate  system.  Thus, 
hda  =  (xfc+i  —  x*)  x  £  di/2  so  that 

j(n)  1  Y  f1  [(**+1  -  Xt)  x  £]  r  (xt+i  +  Xt)  +  t(xt+i  -  Xt)]n+1 

n  +  2^7-1  2  ‘  L  2 

!  "  t\ 

=  (n  +  2)2"+»  ~  X*)  X  ^  *  t(x*+»  +  Xfc)  +  *(x*+i  “  **)]"+1  *  (27) 

In  these  expressions  a  subscript  of  N  +  1  is  equivalent  to  the  subscript  1:  that  is,  Xff+i  =  Xj,. 
Now, 

[(x*+i  -  x*)  x  £]  •  [(x*+i  +  xk)  +  t(xk+i  ~  xt)]  =  [(xt+1  -  xt)  x  £]  •  (xt+1  +  xk) 

=  [(xt+1  +  xk)  x  (xk+i  -  xt)]  •  £ 

=  2(x*  x  Xfc+i)  •  £ 

=  2(*tyt+i  -  Xk+iyk)  (2.8) 

and  therefore, 

lW  =  fj  ^(n  +  2)2‘”W)  ij^*' +  **>  +  ,(x‘+1  *  P-») 

The  tensor  dyadic  (x  +  y)n  can  be  expanded,  but  one  must  be  careful  not  to  use  the 
binomial  theorem  which  assumes  commutativity  of  x  and  y  in  the  terms:  for  example,  xy  ^  yx. 
Rather, 

(x  +  y)n  =  £  ^n-m.m(x,y)  (2.10) 

VRsO 

where  Pn,m(x>y)  is  the  sum  of  all  terms  which  are  permutations  of  n  copies  of  x  and  m  copies 
of  y:  for  example, 

P*.i(x,y)  =  xxy  +  xyx  +  yxx  (2.11) 

The  definition  for  Pn>m  is  extended  to  the  case  Pq.o  by  defining  P0,o  =  1  ■ 


Substitution  of  equation  (2.10)  into  equation  (2.9)  yields 


i**1  - 1,  £  £  p- — <*»+■ + »»■*—  -  *»)■”* 

_  [XkVk+l  -  *k+lVk)  y*  Pn-7m,7m{^k+l  +  Xfc,Xfc+i  ~  Xfc) 

(n  +  2)2»  ^  2m  + 1 


(2.12) 


where  [n/2]  denotes  the  largest  integer  not  exceeding  n/2. 

The  expression  in  equation  (2.12),  though  seemingly  complicated,  yields  simple  expressions 
for  the  first  few  moments.  In  particular, 


l(0)  =  \  X^(**V*+i  “  Zk+iVk) 

1  k=i 
1  N 

I(1)  =  g  5Z(x‘Vk+i  -  *k+iVfc)(x*+i  +  xt) 

°  k=l 
1  N 

l(J)  =  7i  “  x*+l9k)[(xt+i  +  *k)*  +  (*k+l  ~  Xfc)V3] 

1  k=l 
1  N 

=  Y\  13  (**¥*+»  “  *k+i»k) [2xf+1  +  Xfc+iXt  +  XfcXfc+1  +  2x\\ 


(2.13) 


(2.14) 


(2.15) 


Alternatively,  the  components  of  the  centroid  and  the  second  moment  of  area  can  be  written 


4X)  =  g  ]£(*kyfc+i  ~  **+iVk)(*k+i  +  xk) 
t  if 


(2.16) 


(2.18) 


1  J \ 

-  a  Z)(xkyk+1  -  *k+lVk)(»k+i  +  Vk)  (2.17) 

°  t=i 
1  N 

4*  =  T^Z)(x*Vk+l  -  xk+lVk)(xt+l  +  xk+l*k  +  *kl  (218) 

k=i 
1  * 

4*  =  57  ]C(x*Vfc+i  ~  *k+iVk)[2*k+xyt+i  +  n+iyt  +  xkyk+i  +  2xty*]  (2.19) 

k=l 
1  * 

4?  =  To  E(x*Vk+1  -  *k+iyk)[y*+1  +  yk+iyk  +  y*]  (2.20) 

**  k=l 

These  formulae  are  quite  general,  correctly  calculating  the  moments  of  polygons  with  arbi¬ 
trary  connectivity.  For  example,  Figure  1  indicates  a  correct  ordering  of  vertices  for  calculating 
the  moments  of  two  disconnected  squares  while  Figure  2  indicates  a  correct  ordering  of  vertices 
to  be  used  to  calculate  the  moments  of  a  square  containing  a  square  hole.  Interior  holes  must 
be  traversed  clockwise. 


Figure  1:  Order  of  vertices  for  disconnected  polygons. 


2.2  Avoiding  Round-off  Errors 


When  using  the  formulae  derived  in  the  previous  section  for  numerical  calculations,  care 
must  be  taken  to  avoid  round-off  error  in  the  term  (nyt+x  -  xt+iV*)  when  the  origin  of  the 
coordinate  system  is  many  mean  polygon  diameters  from  the  centroid  of  the  panel.  In  this 
case  XkVk+i  »  Xk+iVk  and  the  term  (xkVk+i  -  Zfc+iV*)  is  the  difference  of  two  large,  nearly 
equal  numbers.  Given  any  point  Xo  which  is  close  to  the  polygon,  these  round-off  errors  can  be 
avoided  in  two  ways: 

1.  by  using  the  expression 

XkVk+1  -*fc+i V*  =  (xk  -*o)(yt+i  -  yo)  -  (*k+i  -xo)(Vk  ~  Vo)+xo(Vk+i~Vk)  -Vo(xk+i  -Xk) 

(2.21) 

whose  right  hand  side  does  not  contain  the  product  of  two  large  numbers;  or 

2.  by  first  shifting  the  coordinate  origin  to  xo,  calculating  the  moments,  then  using  the 
formulae 

l(°)  =  40)  (2.22) 

jC1)  —  f  jcdxdy  =  f  (x-xo )dxdy  +  f  Xodxdy  =  +  Xol^  (2.23) 

Ja  Ja  J  a 

l(J)  =  f  x*dxdy 
Ja 

=  /  (x-Xo )7dxdy+  I  xxo dxdy+  /  xo xdxdy  -  I  xodxdy 

Ja  Ja  Ja  Ja 

=  I^J)  +  xoI(1)  +  I(1)xo  -  xSl(0)  (2.24) 

where  the  subscript  0  denotes  a  moment  about  Xo- 

Since  the  second  method  requires  fewer  arithmetic  operations,  it  was  used  in  the  code  provided 
in  Appendix  A.  The  first  vertex  xi  was  chosen  for  Xo. 

Timing  comparisons  on  the  DREA  DEC- 20/60  computer  suggest  that  the  penalty  in  run 
time  paid  for  avoiding  round-off  errors  is  between  10%  and  15%  for  polygons  with  small  number 
of  vertices.  If  speed  is  of  the  essence  and  it  is  known  that  the  origin  of  the  coordinate  system 
will  always  lie  near  the  panel,  it  may  be  best  not  to  worry  about  round-off  errors.  On  the  other 
hand,  the  following  example  (run  on  the  DREA  DEC-20/60)  serves  to  illustrate  the  need  in 
general. 

With  N  =  4,  x\  —  (0,0),  xj  =  (1,0),  xi  =  (1,1),  and  X\  =  (0,1),  the  code  given  in 
Appendix  A  and  similar  code  ignoring  round-off  errors  both  returned  the  following  values  for 
the  moments  (correct  to  7  significant  figures): 

j(°>  =  1.0000000  /W  =  0.5000000  /J1}  =  0.5000000 

li2J  =  0.3333333  =  0.2500000  =  0.3333333 


(2.25) 

(2.26) 


However,  when  each  vertex  was  displaced  by  adding  (10s,  10s)  to  it,  the  code  of  Appendix  A 
returned 

/(°)  =  1.000000  iff  =  1.000005  x  10s  iff  =  1.000005  x  1010  (2.27) 

iff  =  iff  =  iff  =  1.000010  x  1010  (2.28) 

which  are  again  correct  to  7  significant  figures.  However,  the  code  ignoring  round-off  errors 
returned 

/<°>  =  128.00000  41}  =  iff  =  8-566699  x  10®  (2.29) 

iff  =  iff  =  iff  =  6.450141  X  10“  (2.30) 

2.3  Comparison  with  the  Method  of  EN967 

As  mentioned  in  the  introduction,  the  development  of  the  above  formulae  for  the  moments 
of  a  polygon  was  spurred  by  potential  flow  calculations  using  panel  methods.  For  some  years 
DREA  has  used  the  program  EN967  to  calculate  potential  flows.  The  expressions  given  above 
for  the  panel  moments  improve  upon  the  EN967  expressions  in  the  following  ways. 

1.  The  code  is  more  efficient  than  that  used  by  EN967.  Timing  comparisons  indicate  that  the 
new  method  calculates  the  first  three  panel  moments  (the  ones  required  by  the  potential 
flow  calculations)  in  approximately  50%  of  the  time  taken  by  EN967. 

2.  The  new  method  is  more  robust  than  that  used  by  EN967.  The  analytic  expressions 
derived  above,  and  hence  the  code  derived  from  them,  are  completely  free  of  singularities. 
The  EN967  code  relied  upon  manipulations  of  the  slopes  of  the  panel  sides  to  calculate 
the  moments  of  the  panel.  When  the  slopes  were  very  large  or  very  small  (but  non-zero), 
very  large  relative  errors  could  occur  in  the  values  of  the  moments.  As  an  example, 
both  methods  were  used  to  calculate  the  moments  of  the  panel  having  corner  points 
(x,y)  =  (1.0,0.0),(— 1.0,0.5),(— 1.0,— 1.0),  and  (0.0, 0.0).  Both  methods  returned  the 
correct  values  of  the  moments  iff  =  0.41667,  iff  =  0.08333 ,ljff  =  0.10417.  However, 
when  the  third  point  was  changed  to  (—1.0000001,0.5),  the  new  method  returned  the 
correct  values  (which  are  as  before  to  five  significant  figures)  while  EN967  gave  iff  — 
-327,680. 

3.  The  new  method  is  more  general  allowing  arbitrary  polygons.  The  method  of  EN967 
allowed  only  quadrilaterals. 


3  Determining  if  a  Point  is  Inside  a 
Polygon 


The  second  problem  to  be  addressed  is  whether  a  given  point  lies  inside  or  outside  an 
arbitrary  polygon.  It  is  sufficient  to  consider  the  given  point  to  be  the  origin,  since  it  can 
always  be  made  to  be  so  by  a  simple  coordinate  translation. 

Problem  2:  If  x*,  k  =  1,. . .  ,N  are  the  vertices  of  a  polygon  in  order  as  one  proceeds  around 
its  perimeter  counterclockwise,  is  the  origin  inside  or  outside  the  polygon  f 

A  simple  and  efficient  solution  to  Problem  2  can  be  obtained  by  using  the  divergence 
theorem  again.  First  define  the  function  sgn(x)  by 

sgn(x)  =  1  if  x  >  0 

=  0  if  x  =  0  (3.1) 

=  —1  if  x  <  0 


It  has  the  properties 

sgn(x)  =  -sgn(-x) 
sgn(xy)  =  sgn(x)sgn(y) 

sgn  =  sgn(x)  if  x  ^  0 

^sgn(x)  =  2  «(x) 

where  £(x)  is  the  Dirac  delta  function.  It  is  straightforward  to  show  that 

jf  $(*)/{*)<**  =  /(0)[»gn(6)  -  sgn(o)]/2 

and  hence  that 

I  S(x)6(y)dxdy  =1  if  (0,0)  is  inside  the  polygon 

Ja 

=  1/2  if  (0,0)  is  on  an  edge  of  the  polygon 
=  1/4  if  (0,0)  is  on  a  vertex  of  the  polygon 
=  0  if  (0,0)  is  outside  the  polygon 


(3.2) 

(3.3) 

(3.4) 

(3.5) 


(3.6) 


(3.7) 


Using  aquation  (3.5)  and  tba  divergence  theorem  ona  has 

/  f(*)*(v)4*dv  =  f  V  •  (*(*)s«n(y)$jdx<fy 
Jo  Jo 

=  \  faa  M(*)h»(»)^ 

Thai  tima  tha  4th  aida  will  ba  parameterized  with  raapact  to  tba  z  coordinate: 


(*-**) 

\**+i  ~  *kJ 


Than,  n^ds  =  -dx  and,  using  tha  propartiaa  listed  in  aquations  (3.2)-(3.*j,  one  has 

/  *(*)%)***¥ 

J  D 

i£+l  £(x)sgn  [y*  +  ~  **)]  **  i{  **  *  **+» 

0  if  *e  =  z*+i 

(•*n(*s+i)  -  Kn(*k)]sgn  [y*  -  *s]  if  *s  #  *s+i 


'ill 
-ill 
‘ill 

-ill 

1  * 

“7  £|agn(ze+i)  ~  ■*n(z»)]sgn(z*+1  -  i*)«fn(i*+iy*  -  z*y*+i) 
t-i 

1  * 

7  5Z  l>fn(**+i)  -  agn(z*)|agn(z*y*+i  -  *e+iye) 

•  a _ a 


(3.8) 


(3.9) 


0  if  z*  =  z*+j 

(afn(zk+i)  -  sgn(z*)]sgn  if  *s  #  **+i 

0  if  **  =  **+ 

[sgn(ze+i)  -  sgn(x*)]sfn(z*+i  -  **)sfn(z*+1y4  -  zsy*+i)  if  **  *  z*+i 
0  if  z*  =  z*+j 


(310) 


Tha  expression  of  aquation  (3.10)  is  vary  simple  to  calculate.  Note,  however,  that  its 
efficiency  can  ba  enhanced  by  avoiding  unnecessary  multiplications  and  additions;  this  is  done 
in  computer  coda  by  using  appropriate  IF  THEN  ELSE  blocks  or  CASE  statements.  Thus,  an 
efficient  way  to  calculate 


M 

Y.  Ik»(**+i)  -  8gn(z*)|sgn(zey*+i  -  zfc+,y*) 

*»i 


is  by  the  following  algorithm. 


If  Sfc+i  >  0  tkn 

«•  :■  ««■  ♦  3«I*(**lfr+i  -  zft+ift) 
•1m  if  zs+i  =  0  then 

MB  MB  ♦  «C»(«li *+i  -  Zft+ift) 

•*4  if 

•1m  if  si  >  0  tin 
if  Sfc+i  <  0  tin 

MB  :•  MB  ♦  2sf«(x»|f|fc+i  -  **+!»*) 
•1m  if  t^|  =  0  tin 

mb  :•  mb  ♦  ■fa(**V*+i  -  *l+i»i) 
rad  if 

•1m  if  Sft+|  ^  0  tin 

mb  :•  raa  ♦ 

•ad  if 
•ad  do 


Timing  comparisons  oa  the  DREA  DEC-20/60  indicate  that  this  algorithm  is  about  30- 
50%  more  officiant  than  using  aquation  (3.10)  directly.  This  algorithm  is  used  in  the  FORTRAN 
77  subroutine  INPOLY  given  in  Appendix  B. 


Round-off  «rrors  asaociatad  with  tha  term  (****♦ i  ~  **-*■!**)  m  not  as  critical  in  this 
problem  as  in  Problem  1.  In  tha  critical  cam  whan  both  z’s  and  y’s  arc  large,  the  term 
|sgn(z*+i)  -  sgn(xs)|  is  saro,  so  that  there  is  no  contribution  to  tha  sum.  However,  round-off 
errors  could  ba  important  for  points  lying  vary  dam  to  an  edge,  shifting  them  just  enough  so 
that  they  no  longer  lie  insida  (or  outsida)  tha  polygon.  There  is  no  simple  means  of  correcting 
for  this,  but  a  passible  solution  is  to  determine  the  minimum  distance  of  the  origin  to  the 
perimeter,  thus  allowing  the  user  to  decide  whan  tha  origin  is  too  dam  to  an  edge  The 
distance  of  tha  origin  to  tha  4th  aids  is 


!**»*+ i  -  *s+iy*| 

(*i+i  -  **)*  +  (w+i  -  vs) 


(3  11) 


Tha  distance  to  tha  perimeter  is  not  calculated  by  the  subroutine  INPOLY  in  Appendix  B  be¬ 
cause  it  decreases  the  efficiency  of  the  subroutine,  and  is  not  necessary  for  all  uses.  Modification 
of  INPOLY  to  calculate  the  distance  to  the  perimeter  is  straightforward. 


4  Concluding  Remarks 


Analytic  wprwio—  for  the  momenta  of  an  arbitrary  polygon  have  been  derived  and  uaed 
to  develop  a  FORTRAN  77  subroutine  which  calculatee  the  first  three  momenta  of  an  arbitrary 
polygon.  A  similar  expression  (with  corresponding  subroutine)  for  determining  whether  a  point 
lies  within  a  polygon  has  also  been  derived.  These  expressions  have  the  following  properties. 

1.  They  are  simple  and  concise  and  therefore  easily  programmed  in  a  computer  language 
(see  Appendices  A  and  B). 

2.  They  are  computationally  efficient. 

S.  The  code  derived  from  the  expressions  is  robust  provided  care  is  taken  to  avoid  round-off 
errors  in  the  terms  (**y*+i  -  **+iw) 

4.  They  are  general,  providing  correct  expressions  for  polygons  with  any  number  of  sides  or 
any  degree  of  connectivity. 


o  o 


Appendix  A  FORTRAN  77  Subroutine 
PLYMOM 


The  FORTRAN  77  subroutine  PLYMOM  calculates  the  first  three  moments  of  an  arbitrary 
polygon  using  the  methods  discussed  in  Sections  2.1  and  2.2. 

SUBROUTINE  PLYMOM (M. VERTEX. AREA. ACERT.SECMON) 

PLTMQM  calculates  the  first  three  noneate  of  an  arbitrary  polygon. 

C  It  assoaes  counter-clockwise  panel  corner  point  order. 

C 

C  Subroutine  PLYMOM  was  developed  by  the  Canadian  Departaent  of 
C  Motional  Defence. 

C 

C  Author:  David  Bally.  14/1/87 
C 

C  IMPUT: 

C 

CM 

C  VERTEX  • 

C 
C 
C 

C  OUTPUT: 

C 

C  AREA  - 
C  A CERT  • 

C  8ECM0M  - 

C 
C 

INTEGER  M.  MP1.  M 

REAL  AREA.  ACEMT(2) .  SECMOM(S) ,  VERTEX(2.N).  TKXK ,  VSUMI .  VSUMY, 

♦  XKMO.  XKP1M0 .  TXMO.  TKP1N0 

C  Calculate  noneate  about  [VERTEX (1,1), VERTEX (2,1)] 

AREA -00 
ACEMT(l)-0.0 
ACEMT(2)-0  0 


The  nunber  of  vertlcee  of  the  polygon. 

A  2  x  M  array  containing  the  vertices  of  the  polygon. 
VERTEX(l.I)  is  the  x-conponent  of  the  Ith  vertex. 
VERTEK2.I)  is  the  y-conponeat  of  the  Ith  vertex. 


Area  of  the  polygon 

First  nonents  of  the  panel  (centroid  tines  area) 
2nd  noneate  of  area  of  the  panel.  8ECM0N( 1 )-Ixx . 
SECM0N(2)-Ixy-Iyx.  8ECM0N(3)-Iyy 


«tCMON(l)-0.0 

ncM0M(a)-o.o 

•1001(3) -0.0 
ZKNO-O.O 
TKNO-O.O 
DO  10  N-l.l 
MP1-N+1 

IF  (N.EQ.I)  NF1-1 
XXP1M0-VERTEX(1  ,NP1)  -  VERTEX  (1 . 1) 

TKP1M0-VERTEX (2 . MP1 ) - VERTEX (2 . 1 ) 

VSUMX-XKP1NO+XKMO 

VSUNY-YKP1N0-YKM0 

C  Calculate  AREA 

TKIK-TKP1MO*XKNO'TKMO*XKP1MO 

AREA-AREA+TKXK 

C  Calculate  ACEHT 

ACEHT ( 1 ) -ACEMT( 1 ) ♦VSOMX*YKXK 
ACEHT (2) -ACEHT(2) +VSUMY-YKXK 

C  Calculate  SECMOM 

SECNON ( 1 ) -SECMOM (1 ) +TKXK* (XKP 1M0-VSUMX+XKM0*  *2) 

SECMOM (2) -SECMOM (2) +YKIK* (XKP1M0*YKP1M0+XKM0*YKM0+VSUMX*VSUMY) 
SECMQM(3)-SECM0M(3)+TKXK*(YKP1M0*VSTJMY+YKM0**2) 

XXM0-XKP1M0 
TKM0-TKP1M0 
10  COITIIUE 

C  Calculate  aoueats  about  (0,0) 

AREA-AREA-0.5 

ACERT ( 1 ) -ACEHT ( 1 ) /6 . 0- VERTEX ( 1 , 1 )- AREA 
ACEHT (2) -ACEHT (2) /O . 0- VERTEX (2,1) -AREA 
SECNON ( 1 ) -8ECM04  ( 1 )/ 1 2 . 0- VERTEX ( 1 , 1 )*( 2 . 0-ACENT ( 1 ) - 

-  VERTEX (1.1) -AREA) 

8ECN0M(2)-8ECM0M(2)/24.O-VERTEX(l ,1)*ACENT(2)+ 

-  VERTEX(2, 1)*(ACEHT(1) ‘VERTEX (1 , 1)-AREA) 

SECMOM (3) -SECMOM (3) /12 . 0-VERTEX(2 , 1) * (2 .0-ACENT(2) - 

♦  VERTEX (2.1) -AREA) 

RETDRH 

EHD 


Appendix  B  FORTRAN  77  Subroutine 
INPOLY 


The  FORTRAN  77  subroutine  INPOLY  (unction  determines  whether  a  point  is  inside  an 
arbitrary  polygon  using  the  method  discussed  in  Section  S. 

SUBROUTINE  INPOLY (N .VERTEX , X , IFLAG) 

C 

C  INPOLY  determines  whether  the  point  X  lies  inside  an  arbitrary 
C  polygon.  It  assuaes  counter-clockwise  panel  corner  point  order. 

C 

C  Subroutine  INPOLY  was  developed  by  the  Canadian  Department  of 
C  National  Defence. 

C 

C  Author:  David  Hally.  14/1/87 
C 

C  INPUT: 

C 

C  N  *  The  nunber  of  vertices  of  the  polygon. 

C  VERTEX  •  A  2  x  N  array  containing  the  vertices  of  the  polygon. 

C  VERTEX (1. I)  is  the  x-coaponent  of  the  Ith  vertex. 

C  VERTEX(2,I)  is  the  y-coaponent  of  the  Ith  vertex. 

C  X  -An  array  of  length  2  containing  the  point  which  is  to  be 

C  checked. 

C  X(l)  is  the  x-coaponent  of  the  point. 

C  X(l)  is  the  y-coaponent  of  the  p~int. 

C 

C  OUTPUT: 

C 

C  IFLAG  *  4,  if  X  is  Inside  the  polygon 

C  ■  2,  if  X  is  on  s  side  of  the  polygon 

C  •  1,  if  X  is  on  s  vertex  of  the  polygon 

C  ■  0,  if  X  is  outside  the  polygon 

C 

INTEGER  IFLAG.  N.  MP1 .  N.  SGN 

REAL  VERTEX (2, N).  X(2) ,  XKMO.  XKP1N0,  YKMO,  YKP1M0 
I FLAG-0 


14 


XKM0-VERTEX(1.1)-X(1) 

YKHO-VERTEX(2 . 1) -X(2) 

DO  10  N-l.N 
MP1-M+1 

IF  (N.EQ.H)  MP1-1 

*  XKP1M0-VERTEX(1.MP1)-X(1) 

TKP1M0-VERTEX(2 ,MP1) -X(2) 

IF  (XKMO.LT.O)  THEN 

-  IF  (XKP1M0.GT.0)  THEN 

IFLAG-IFLAG+2*SGN(XKM0*YKP1M0-XKP1M0*YKM0) 
ELSE  IF  (XXP1M0.EQ.0)  THEN 

IFLAG-IFLAG+SGN (IKM0*YKP1M0-IKP1M0*YKM0) 
END  IF 

ELSE  IF  (XKMO.GT.O)  THEN 
IF  (XKP1M0.LT. 0)  THEN 

IFLAG-IFLAG+2*SGN(XKM0*YKP1M0-XKP1M0*YKM0) 
ELSE  IF  (XKP1M0.EQ .0)  THEN 

IFLAG-IFLAG+SGN (IKMO+TKP1MO-IKP1MO+TKNO) 
END  IF 

ELSE  IF  (XKP1M0.NE.0)  THEN 

IFLAG-IFLAG+SGN (XKMO+YKP 1M0-XKP 1M0+TKM0) 

END  IF 

XKM0-XKP1M0 

TKM0-TKP1N0 


10  CONTINUE  - 

RETURN  \ 

END 

INTEGER  FUNCTION  SGN(X) 

C  Calculates  sgn(X) 

REAL  X 

IF  (X.GT.O.O)  THEN 
8GN-1 

ELSE  IF  (X.LT.O.O)  THEN 

8GN— 1  | 

ELSE 

8GN-0  i; 

END  IF 
RETURN 
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